#delimit;
clear all;
version 13.1;
pause on;
program drop _all;
capture log close;
set more off;

/*INSERT FOLDER PATH*/
cd /;


use pmra_level_freq.dta, clear;
gen yr_to_rtrct=year-rtrct_year;

forvalues i=5(-1)1 {;
	gen yrs_before_rtrct`i'=0;
	replace yrs_before_rtrct`i'=1 if yr_to_rtrct==-`i' & rtrct==1;
};


forvalues i=1(1)10 {;
	gen yrs_after_rtrct`i'=0;
	replace yrs_after_rtrct`i'=1 if yr_to_rtrct==`i' & rtrct==1;
};
gen yrs_after_rtrct11=0;
replace yrs_after_rtrct11=1 if yr_to_rtrct>10 & rtrct==1;

gen yrs_before_rtrct6=0;
replace yrs_before_rtrct6=1 if yr_to_rtrct<-5 & rtrct==1;

drop yr_to_rtrct;

keep id year nb_pmra rtrct rtrct_year treat case_code yrs_before_rtrct6 yrs_before_rtrct5-yrs_after_rtrct11;
order id year nb_pmra rtrct rtrct_year treat case_code yrs_before_rtrct6 yrs_before_rtrct5-yrs_after_rtrct11;

xtqmlp nb_pmra yrs_before_rtrct6-yrs_before_rtrct1 yrs_after_rtrct1-yrs_after_rtrct11 i.year, fe i(id) cluster(case_code);
estimates save estimates/poisson_freq_dynamics.ster, replace;


estimates drop _all;
estimates use estimates/poisson_freq_dynamics.ster;

matrix B=e(b);
matrix V=e(V);

keep if rtrct==1;
gen t_to_rtrct=year-rtrct_year;
bysort t_to_rtrct: keep if _n==_N;

foreach y in 1 2 3 4 5 {;
	local t=colnumb(B,"yrs_before_rtrct`y'");
	gen yrs_before_rtrct`y'_coef = B[1,`t'];
	gen yrs_before_rtrct`y'_var = V[`t',`t'];
};

foreach y in 1 2 3 4 5 6 7 8 9 10 {;
	local t=colnumb(B,"yrs_after_rtrct`y'");
	gen yrs_after_rtrct`y'_coef = B[1,`t'];
	gen yrs_after_rtrct`y'_var = V[`t',`t'];
};

gen treatment_coefs=0;
gen treatment_vars=0;

foreach var in 	yrs_before_rtrct5 yrs_before_rtrct4 yrs_before_rtrct3 yrs_before_rtrct2 yrs_before_rtrct1 
 		yrs_after_rtrct1 yrs_after_rtrct2 yrs_after_rtrct3 yrs_after_rtrct4 yrs_after_rtrct5 
 		yrs_after_rtrct6 yrs_after_rtrct7 yrs_after_rtrct8 yrs_after_rtrct9 yrs_after_rtrct10 
 		{;
		replace treatment_coefs=`var'_coef  if `var'==1;
		replace treatment_vars=`var'_var  if `var'==1;
		};

gen treatment_95low=treatment_coefs-1.96*sqrt(treatment_vars);
gen treatment_95hi=treatment_coefs+1.96*sqrt(treatment_vars);

sort t_to_rtrct;

twoway scatteri -2.00 0 2.00 0, c(l) lcolor(black) lwidth(thin) msym(none) yline(0, lcolor(black) lwidth(thin))
|| line treatment_95hi t_to_rtrct, lpattern(shortdash) lcolor(red) lwidth(medium) || line treatment_95low t_to_rtrct, 
lpattern(shortdash) lcolor(red) lwidth(medium) || line treatment_coefs t_to_rtrct, lpattern(solid) lcolor(blue) lwidth(medthick) || 
if t_to_rtrct >= -5 & t_to_rtrct <= 10, xtitle(" " "Time to/After Retraction") xlabel(-5(1)10, labsize(small)) ylabel(-2.00(0.50)2.00, 
format(%15.2f) angle(horizontal) labsize(small)) legend(off) saving(graphs/poisson_freq_dynamic.gph, replace);




/*funding*/

use pmra_level_funding.dta, clear;
replace amount=amount/1000000;
gen yr_to_rtrct=year-rtrct_year;

forvalues i=5(-1)1 {;
	gen yrs_before_rtrct`i'=0;
	replace yrs_before_rtrct`i'=1 if yr_to_rtrct==-`i' & rtrct==1;
};


forvalues i=1(1)10 {;
	gen yrs_after_rtrct`i'=0;
	replace yrs_after_rtrct`i'=1 if yr_to_rtrct==`i' & rtrct==1;
};
gen yrs_after_rtrct11=0;
replace yrs_after_rtrct11=1 if yr_to_rtrct>10 & rtrct==1;

gen yrs_before_rtrct6=0;
replace yrs_before_rtrct6=1 if yr_to_rtrct<-5 & rtrct==1;

drop yr_to_rtrct;

keep id year amount rtrct rtrct_year treat case_code yrs_before_rtrct6 yrs_before_rtrct5-yrs_after_rtrct11;
order id year amount rtrct rtrct_year treat case_code yrs_before_rtrct6 yrs_before_rtrct5-yrs_after_rtrct11;

xtqmlp amount yrs_before_rtrct6-yrs_before_rtrct1 yrs_after_rtrct1-yrs_after_rtrct11 i.year, fe i(id) cluster(case_code);
estimates save estimates/poisson_fndng_dynamics.ster, replace;


estimates drop _all;
estimates use estimates/poisson_fndng_dynamics.ster;

matrix B=e(b);
matrix V=e(V);

keep if rtrct==1;
gen t_to_rtrct=year-rtrct_year;
bysort t_to_rtrct: keep if _n==_N;

foreach y in 1 2 3 4 5 {;
	local t=colnumb(B,"yrs_before_rtrct`y'");
	gen yrs_before_rtrct`y'_coef = B[1,`t'];
	gen yrs_before_rtrct`y'_var = V[`t',`t'];
};

foreach y in 1 2 3 4 5 6 7 8 9 10 {;
	local t=colnumb(B,"yrs_after_rtrct`y'");
	gen yrs_after_rtrct`y'_coef = B[1,`t'];
	gen yrs_after_rtrct`y'_var = V[`t',`t'];
};

gen treatment_coefs=0;
gen treatment_vars=0;

foreach var in 	yrs_before_rtrct5 yrs_before_rtrct4 yrs_before_rtrct3 yrs_before_rtrct2 yrs_before_rtrct1 
 		yrs_after_rtrct1 yrs_after_rtrct2 yrs_after_rtrct3 yrs_after_rtrct4 yrs_after_rtrct5 
 		yrs_after_rtrct6 yrs_after_rtrct7 yrs_after_rtrct8 yrs_after_rtrct9 yrs_after_rtrct10 
 		{;
		replace treatment_coefs=`var'_coef  if `var'==1;
		replace treatment_vars=`var'_var  if `var'==1;
		};

gen treatment_95low=treatment_coefs-1.96*sqrt(treatment_vars);
gen treatment_95hi=treatment_coefs+1.96*sqrt(treatment_vars);

sort t_to_rtrct;

twoway scatteri -2.50 0 2.50 0, c(l) lcolor(black) lwidth(thin) msym(none) yline(0, lcolor(black) lwidth(thin))
|| line treatment_95hi t_to_rtrct, lpattern(shortdash) lcolor(red) lwidth(medium) || line treatment_95low t_to_rtrct, 
lpattern(shortdash) lcolor(red) lwidth(medium) || line treatment_coefs t_to_rtrct, lpattern(solid) lcolor(blue) lwidth(medthick) || 
if t_to_rtrct >= -5 & t_to_rtrct <= 10, xtitle(" " "Time to/After Retraction") xlabel(-5(1)10, labsize(small)) ylabel(-2.50(0.50)2.50, 
format(%15.2f) angle(horizontal) labsize(small)) legend(off) saving(graphs/poisson_fndng_dynamic.gph, replace);